function wt_coef_gal = wavelet_OMLSA_cohen(wt_coef,lambdaD,alpha)
% This function is exactly by Cohen's 2001 wavelet paper
size_wt=size(wt_coef);
num_ite=max(size_wt(1),size_wt(2));
beta=0.7;
qmax=0.95;
ksimin=0.1;
ksimax=10^(-0.5);

for I=1:num_ite
  gamaK = wt_coef(I)^2/lambdaD(I);
   if I==1
      ksi(I)=alpha+(1-alpha)*p(gamaK-1);
      TEMP=ksi(I)/(ksi(I)+1);
      wt_coef_gal(I)=TEMP*exp(0.5*expint(TEMP*gamaK))*abs(wt_coef(I));
      ita(I)=ksi(I);
      ksimin=ksi(I);
      ksimax=ksi(I)+0.0001;
      
   else
      ksi(I)=beta*ksi(I-1)+(1-beta)*ita(I-1);
      temp_local=log10(ksi(I)/ksimin)/log10(ksimax/ksimin);

      P_local=(ksi(I)>=ksimax)+temp_local.*...
                      (ksimin<ksi(I)&ksi(I)<ksimax);
      q=1-P_local;
      q = qmax*(q>=qmax)+q.*(q<qmax);
      %========== equation (10) ================
      ita(I)=alpha*abs(wt_coef_gal(I-1))+(1-alpha)*p(abs(wt_coef(I))-sqrt(lambdaD(I)));  
      lambdaX=ita(I)*lambdaD(I);
    
      %=========================================
      v=gamaK.*ita(I)./(1+ita(I));
      % The only difference with Cohen's paper is without 'sqrt'
      pp=1./(1+q*sqrt(1+ita(I)).*exp(-v./2)./(1-q));
      %pp=1./(1+q*(1+ita(I)).*exp(-v./2)./(1-q));
      wt_coef_gal(I)=(lambdaX.*pp)./(lambdaX+lambdaD(I)).*wt_coef(I);      
      
   end
end

return
